clear all
close all

% =====================================================================
% ====  "Risk-Sharing in Village Economies Revisited", by Bold and Broer

% ====   This programme solves the limited commitment village economy for
% both the standard specification (individual deviations) and the
% alternative specification (coalitional deviations)

% ====    Written by Tessa Bold and Tobias Broer

% ====    This version June 2021
% =====================================================================


filetosave='Vilcoal_03_02_21_vill_cd_incproc_noequal'; % name of workspace saved

%====================================================
% select parameters
% ===================================================

village=1;                         % village indicator: 1 - Aurepalle; 2-Kanzara; 3-Shirapur
coaldev=1;                         % set to 1 for coalitional deviations, to 2 for self-insurance (SI), to 0 for individual deviations (ID), to 3 for full-insurance with coalitional deviations
Tessa1Tobi2=1;                     % sets paths depending on user below
unconditional=0;                   % set this to 1 for unconditional income data, to 0 for income data conditional on time fixed effets / demeaned data period-by-period
fixedeffect=1;                     % set this to 1 for an income process with fixed effects
one_per_aut=0;                     % set this to 1 if you want deviations to imply one period of autarky
ytaxrate=0;                        % linear tax on income
scatter = 0;		           % set to 1 to just draw a scatter plot
Rov_sb = 0;                          % set this to 1 if you want the outside option for the Rov in the case of coalitional deviations to equal the average value from the solution to the next smaller village (second-best), rather than
                                   % just the utility from consuming average ROV income (first-best)
estrho1=0;			   % set this to 1 if you would like to estimate persistence, default is 0
rhogrid=0;                         % set to 1 if you want to have also a grid for CRRA rho - only to show non-identification in Aurepalle (village=1)
Nphi=1;                            % number of borrowing limits in SI economy - not used currently
server=0;                          % set to 1 to run on server; o wise path for workstations is used
buildup=0;                         % set to 1  to calculate ID model for subset of group sizes, not just whole village
Nincproc=1;                       % number of income processes to look at
if estrho1==0
maxincerr=2/3;                     % maximum income measurment error in multiples of income variance (in levels)
elseif estrho1==1
maxincerr=0;
end
model_var_dy=1;                    % set to 1 if you want to use the model-implied variance of dlog(y) to specify maximum meas error, otherwise var(dlog(y)) from data is used, slightly different as we target var(log(y)), not var(dlog(y())
T=6200; Tinit=201;                 % number of simulation periods in total, T-Tinit+1 equals number of initial periods to discard
maxcerr=0.9;                       % max cons meas error in log levels
if scatter==0
N_cerr=0;                         % number of support points for cons meas error
elseif scatter==1
N_cerr=0;
end
if coaldev==0
    convlength=0;                  % the conv criterion must be fulfilled at least in convlength consecutive periods to avoid jumping 'in'
else
    convlength=15;
end
if village==1
    vss_high_inddev=34-1;          % the number of total village members is vss_high_inddev+1 under individual deviations
    vss_high_coaldev=34-1;
elseif village==2;
    vss_high_inddev=37-1;                 
    vss_high_coaldev=37-1;
elseif village==3
    vss_high_inddev=31-1;                
    vss_high_coaldev=31-1;
end
vssbuildupvec_inddev=[1:7,10,15,vss_high_inddev]; % for buildup==1, this is the vector of village sizes for which the solution is calculated in the individual deviations case
if server==1
    if unconditional==0
        outputpath=sprintf('/home/tbroer/Tessa/Output estimation May 2019 conditional');
    elseif unconditional==1
        outputpath=sprintf('/home/tbroer/Tessa/Output estimation May 2019 unconditional');
    end
    programpath=sprintf('/home/tbroer/Tessa');
else
    if Tessa1Tobi2==2
        outputpath=sprintf('C:/Users/tbroer/Dropbox/Research/Current_Projects/tessa/Programmes/Output estimation August 2018 conditional')
        programpath=sprintf('C:/Users/tbroer/Dropbox/Research/Current_Projects/tessa/Programmes')
    elseif Tessa1Tobi2==1
        outputpath=sprintf('C:/Users/tbold/Dropbox/Tessa and Tobi/JEEA Replication/Programmes')
        programpath=sprintf('C:/Users/tbold/Dropbox/Tessa and Tobi/JEEA Replication/Programmes')
    end
end
        cd(programpath)
if buildup==0
    if rhogrid==0 && scatter == 0
        betavec0=[0.5:0.04:0.78,0.8:0.01:0.98];
elseif scatter ==1
betavec0=[0.83];
    elseif rhogrid==1
        betavec0=[0.83,0.85,0.87,0.89,0.91,0.93,0.95,0.97];
    end
elseif buildup==1

betavec0=[0.86,0.905,0.91,0.94]; %betavec0=[0.7:0.04:0.82,0.84:0.02:0.98];

end
if rhogrid==0
if scatter==1
betavec1=[0.94];
else
 betavec1=[0.84:0.02:0.88,0.89:0.01:0.96,0.9625:0.0025:0.99];
 betavec1=[0.95:0.01:0.98];
end
elseif rhogrid==1
    betavec1=[0.84,0.86,0.88,0.90,0.92,0.94,0.96,0.98];
end
betavec3=[0.96:0.0025:0.99];
if coaldev<2 || coaldev==3
    
    rhovec=[1];
elseif coaldev==2
    rhovec=ones(1,Nphi);
end

if rhogrid==0
    gapcritset=0.0001;              %convergence criterion for change in values
elseif rhogrid==1
    gapcritset=0.001;
end

rov_average=1;                      % set to 1 to specify the rest of village in terms of average
                                    % consumption
momentclevinlog=1;                  % set this to 1 if you want to calculate the moments for levels of
                                    % consumption based on log consumption
maxitgapset=400;                    % maximum number of iterations
Pun_fac=0.0;                        % consumption penalty in the first period after deviation equal
                                    % to caut*Pun_fac
% grid of time-varying planner weights
lambda=[0.05:0.01:0.95]';
m=length(lambda);
if floor((m+1)/2)~=((m+1)/2)
    disp('care: lambda has to have uneven no of points in first dimension')
    stop
end
S_KEEP=nan(1000*3,2,vss_high_inddev);               % 
S_rov_KEEP=nan(100,1,vss_high_inddev);              % this one is for a maximum of vss=19 for the CD
P_rov_KEEP=nan(100,1000,vss_high_inddev);
PHI_KEEP=nan(m,1000*3,2,vss_high_inddev);
IDIO_DIST_ROV_KEEP=nan(1000*3,3,vss_high_inddev);

% ====================================================
% Specify Income process for individual and Rest of Village (rov)
% ====================================================
% data moments for log income from ICRISAT data - all villages
% first col: Aurepalle, second: Kanzara, third: Shirapur
if fixedeffect==1
    if unconditional==0
    datafile='armoments_fe0';
elseif unconditional==1
    datafile='armoments_fe1';
    end
elseif fixedeffect==0
     datafile='armoments';
end
if server==0 & Tessa1Tobi2==2
    cd('C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\JEEA Replication\Data\Output')
elseif server==0 & Tessa1Tobi2==1
    cd('C:\Users\tbold\Dropbox\Tessa and Tobi\JEEA Replication\Data\Output')
end
eval(['load ' datafile '.out;']);
eval(['covyylag=' datafile '(1,:);']);
eval(['vardy=' datafile '(2,:);']);
eval(['vary=' datafile '(3,:);']);
if model_var_dy==1
    vardy=[0.231,0.125, 0.269];
end
cd(programpath)
if maxincerr>0
    if Nincproc>1
        varnu_max=maxincerr*vardy/2;                % nu is measurement error in log levels
    else
        varnu_max=zeros(1,3);
    end
    rho_min=covyylag(village)/vary(village);                      % AR(1) coefficient w/o meas error
    rho_max=covyylag(village)/(vary(village)-varnu_max(village)); % AR(1) coefficient w max meas error
    rho1vec=linspace(rho_min,rho_max,Nincproc);
    varnu=max(vary(village)-covyylag(village)./rho1vec,0); % max prevents the first entry to be -e17
    vareps=(vary(village)-varnu).*(1-rho1vec.^2);
else
    rho1_data=covyylag(village)/vary(village);
    vareps_data=vary(village).*(1-rho1_data.^2);
    rho1vec=[-0.6,-0.4,-0.2,0,0.2,0.4,0.6];
    vareps=vary(village).*(1-rho1vec.^2);
    varnu=zeros(1,length(vareps));              
    Nincproc=length(vareps);
end

% ===================================
% Iterate over income processes
% ===================================
    disp('now solving for village,coaldev, unconditional, Nincproc, fixedeffect one_per_aut Rov_sb - press any key to continue')
    disp([village coaldev unconditional Nincproc fixedeffect one_per_aut Rov_sb])
    pause
indicator1=zeros(1,size(rhovec,2));
incproc=1
    
 filetosaveest=['Lambda'];

  
    % specify a markov process for individual income
    [incomevals, Q1] = rouwenhorst(rho1vec(incproc),vareps(incproc)^0.5,3);   
    incomevals=exp(incomevals');
    QQQ=Q1^1000;
    SS=QQQ(1,:)';
    incomevals=incomevals/(SS'*incomevals);
    cumQ=cumsum(Q1,2);
    incomevals=(1-ytaxrate)*incomevals+ytaxrate; % tax as in Krueger and Perri JET 2011
    
    % set minimum and maximum village size
    if coaldev==1                       
        if rhogrid==0
            vss_set=[1:7];% prior solutions have shown that 1-7 is the relevant set. Otherwise use, for example, [1:12,17,22,27,vss_high_coaldev]
        elseif rhogrid==1
            vss_set=[1:5]
        end
        vss_high=vss_set(end);
    elseif coaldev==0
        if buildup==0
            vss_set=[vss_high_inddev];
            lagind=1;
            vss_high=vss_set(end);
        elseif buildup==1
            vss_set=vssbuildupvec_inddev;
            lagind=1;
            vss_high=vss_set(end);
        end
    elseif coaldev==2;
        vss_set=vss_high_inddev;
        vss_high=vss_set(end);
    elseif coaldev==3
        vss_set=[1:10];
        vss_high=10;
    end
    
    
    % ===================================
    % iterate over discount factors
    % ===================================
    if coaldev==0
        betavec=betavec0;
    elseif coaldev==1
        betavec=betavec1;
    elseif coaldev==2
        r=0.04;
if scatter==0
        betamax=1/(1+r)-0.015;
        betavec=linspace(0.8,betamax,20);
elseif scatter ==1
betavec=0.92;
end
        Z=incomevals;
        phinat=min(Z)/(r-1);
        if Nphi>1
            phivec=linspace(phinat,0,Nphi);
        else
            phivec=0;
        end
    elseif  coaldev==3
        betavec=betavec3;
    end
    Momentsest=nan(35,N_cerr+1,size(betavec,2),1,vss_high_inddev,2);
    for equal=0:1
    for Qbeta=1:size(betavec,2)
        beta=betavec(Qbeta);
        if rhogrid==1
            if coaldev==0
                if beta==0.83
                    rhovec=1.468461538461538
                elseif beta==0.85
                    rhovec=1.316666666666667
                elseif beta==0.87
                    rhovec=1.16
                elseif beta==0.89
                    rhovec=1
                elseif beta==0.91
                    rhovec=0.8340
                elseif beta==0.93
                    rhovec=0.6611
                elseif beta==0.95
                    rhovec=0.4808
                elseif beta==0.97
                    rhovec=0.2912
                end    
            elseif coaldev==1
                if beta == 0.98
                    rhovec=0.366222222222222
                elseif beta==0.84
                    rhovec=2.389
                elseif beta==0.86
                    rhovec=2.11
                elseif beta==0.88
                    rhovec=1.821
                elseif beta==0.90
                    rhovec=1.5515
                elseif beta==0.92
                    rhovec=1.28718
                elseif beta==0.94
                    rhovec=1
                elseif beta==0.96
                    rhovec=0.6935
                else
                end
            end
        end

        % ===================================
        % iterate over CRRA
        % ===================================
        for Qrho=1:size(rhovec,2)
            rho=rhovec(Qrho)  ;
            if coaldev==1
                lagind=[];
            else
                lagind=1;
            end
            
            P_LAG=nan(1000,1000,length(vss_set));Y_lag=nan(1000,1000,length(vss_set));
            waut_lag=[]; waut_avg_lag=[]; waut_rov_sb=[]; waut_rov=[];waut_temp=[]; waut_lag_lambda=[];
            WAUT_LAG=nan(1000,2,length(vss_set));
            WAUT_LAG_LAMBDA=nan(length(lambda),1000,2,length(vss_set));
            STABLE=[];
           if coaldev==3
                Stable=nan(length(vss_set));
                stable=[]; countstable=0;
                clear Delta
           end
            if Qbeta<=6
               vss_set=[1 2 3 4]; 
            else
               vss_set=[1 2 3 4 5 6];
            end
            for vsscount=1:length(vss_set)       % iterate over village size
                
                vss=vss_set(vsscount);
                vss_max=max(vss_set);
                if vsscount>1 & coaldev==1
                    lagind=lagind+(vss-(vss_set(vsscount-1)+1));
                end
                disp('now solving for parameters')
                disp([Qbeta,Qrho])
                disp('vss')
                disp(vss)
               
                    % this is the scaling factor: if you want RoV consumption as average consumption per HH, scale rov variables
                    % by vss
                    if rov_average==1
                        ScaleRov=vss;
                    else
                        ScaleRov=1;
                    end
                    csim=[];  cc =[]; eeewaut =[]; eewaut =[]; eewaut_lambda =[]; ewaut =[]; waut =[]; waut_check=[]; waut_check_ind=[]; rr =[]; cshare=[]; raw =[]; cshare =[]; dcshare =[]; dcshare_raw =[]; S_rov =[]; P_rov =[];
                    S_rov_new =[]; P_rovnew =[]; S =[]; P =[]; caut =[]; c =[]; waut =[]; gammagrid =[]; gammar =[]; phi =[]; phisol =[]; csol =[];csolold=[]; cfirstbest =[];
                    wfirstbest =[]; ewfirstbest =[]; wsbnewsol=[];
                    Y =[]; vs=[];  idio_dist_autnew=[];idio_dist_aut_new=[];normperf=[];w1=[];indSj=[]; indSjj=[]; Ssim=[]; aggvillincome=[];aggincgrowth=[];aggvillincgrowth=[];
                    S_rov_aut=[]; S_rov_sb=[]; S_rov_aut_all=[];S_rov_aut_all2=[];
                    w =[];  phinew =[]; Yagg =[]; wint =[]; wsb =[]; ewsb =[]; csol =[]; phisol =[]; nonconv =[];   gammarfun =[];  x0 =[]; index =[]; ewsb =[]; bind=[];
                    waut=[];
                    
                    % ===================================
                    % build the state space of ROV income
                    % ===================================
                    idio_dist=[];
                    idio_dist_aut=[];
                    idio_dist_aut_all=[];idio_dist_aut_all2=[];
                    idio_dist_sb=[];
                    
                    for k1=0:vss
                        for k2=0:vss
                            for k3=0:vss
                                II=zeros(1,3);
                                if k1+k2+k3==vss
                                    % the vector of aggregate incomes from
                                    % permutations of vss households
                                    S_rov=[S_rov; k1*incomevals(1,1) + k2*incomevals(2,1) + k3*incomevals(3,1)];
                                    % the corresponding distribution of
                                    % individual income values
                                    idio_dist=[idio_dist; k1 k2 k3];
                                    
                                    % this computes the 'best' distribution of incomes
                                    % of size vss-lagind
                                    if vss>1
                                        II=zeros(1,3);
                                        temp=[k1 k2 k3];
                                        
                                        
                                        for jj=1:lagind
                                            II(min(find(temp>0)))=II(min(find(temp>0)))+1;
                                            temp=[k1 k2 k3]-II;
                                        end
                                        idio_dist_aut=[idio_dist_aut; [k1 k2 k3]-II];
                                        S_rov_aut=[S_rov_aut; ([k1 k2 k3]-II)*incomevals];
                                        
                                        
                                        II=zeros(1,3);
                                        temp=[k1 k2 k3];
                                        
                                        for jj=1:lagind-1
                                            II(min(find(temp>0)))=II(min(find(temp>0)))+1;
                                            temp=[k1 k2 k3]-II;
                                        end
                                        idio_dist_sb=[idio_dist_sb; [k1 k2 k3]-II];
                                        S_rov_sb=[S_rov_sb; ([k1 k2 k3]-II)*incomevals];
                                    end
                                end
                            end
                        end
                    end
                    
                    % resort S_rov in ascending order
                    [S_rov,S_rovind]=sort(S_rov);
                    idio_dist=idio_dist(S_rovind,:);
                    
                    if vss>1
                        idio_dist_aut=idio_dist_aut(S_rovind,:);
                        S_rov_aut=S_rov_aut(S_rovind);
                        idio_dist_sb=idio_dist_sb(S_rovind,:);
                        S_rov_sb=S_rov_sb(S_rovind);
                    end
                    
                    
                    %%%%and create S_rov_aut vectors that drops lagind individuals from
                    %%%%the rov vector. Why lagind? The biggest stable group is
                    %%%%vss+1-lagind, hence, there will be vss+1-lagind-1
                    %%%%invididuals from the rov in it, i.e. vss-lagind, hence
                    %%%%lagind individuals will not be in it. 
                    if vss>1
                        for k=1:length(STABLE)
                            if STABLE(k)==1
                                for g=1:length(idio_dist)
                                    index=0;
                                    for j1=0:vss-k
                                        for j2=0:vss-k
                                            for j3=0:vss-k
                                                if j1+j2+j3==vss-k && (idio_dist(g,1)>=j1 & idio_dist(g,2)>=j2 & idio_dist(g,3)>=j3)
                                                    index=index+1;
                                                    temp=[j1 j2 j3];
                                                    idio_dist_aut_all(g,:,index,k)=[idio_dist(g,:)-temp];
                                                    S_rov_aut_all(g,:,index,k)=[(idio_dist(g,:)-temp)*incomevals];
                                                end
                                            end
                                        end
                                    end
                                end
                            end
                            
                        end
                    end
                    
                    % ===================================
                    % simulate aggregate income transitions
                    % ===================================
                    tic
                    T_trans=50000;
                    aggincind=zeros(length(S_rov),length(S_rov));
                    for jjj=1:length(S_rov)
                        jjj;
                        inccount=zeros(T_trans,vss,3);
                        agginc=zeros(T_trans,3);
                        for iii=find(idio_dist(jjj,:)>0)
                            for kkk=1:idio_dist(jjj,iii)
                                rng(jjj*100+iii*10+kkk,'twister');
                                shock=rand(T_trans,1);
                                rr=(find(shock<=cumQ(iii,1)));
                                inccount(rr,kkk,iii)=1;
                                rr=(find(cumQ(iii,1)<shock & cumQ(iii,2)>=shock));
                                inccount(rr,kkk,iii)=2;
                                rr=(find(cumQ(iii,2)<shock & cumQ(iii,3)>=shock));
                                inccount(rr,kkk,iii)=3;
                            end
                            agginc(:,iii)=sum(incomevals(inccount(:,1:idio_dist(jjj,iii),iii)),2);
                        end
                        agginc=sum(agginc,2);
                        for jjj2=1:length(S_rov)
                            qqq=(agginc<=S_rov(jjj2)+0.0000000001).*(agginc>=S_rov(jjj2)-0.0000000001);
                            aggincind(jjj,jjj2)=sum(qqq);
                        end
                        P_rov(jjj,:)=aggincind(jjj,:)/sum(aggincind(jjj,:));
                    end
                    disp('Simulation for agg income process takes')
                    toc

                    
                    if vss>1 && coaldev<2
                        S_rov_aut=S_rov_aut/(max(ScaleRov-lagind,1));
                        S_rov_sb=S_rov_sb/(max(ScaleRov-lagind+1,1));
                        for g=vss_set(find(vss-1>=vss_set))
                            if STABLE(g)==1
                                S_rov_aut_all(:,:,:,g)=S_rov_aut_all(:,:,:,g)/(g);
                            end
                        end
                    end
                    S_rov=S_rov/ScaleRov;
                    
                    % ===================================
                    % state space: combinations of indiv and village income
                    % ===================================
                    S=[kron(incomevals,ones(size(S_rov))),kron(ones(size(incomevals)),S_rov)];
                    idio_dist_rov=kron(ones(size(incomevals)),idio_dist);
                    idio_dist_aut=kron(ones(size(incomevals)),idio_dist_aut);
                    S_rov_aut=[kron(ones(size(incomevals)),S_rov_aut)];
                    S_rov_sb=[kron(ones(size(incomevals)),S_rov_sb)];
                    
                    if vss>1 && coaldev<2
                        for g=vss_set(find(vss-1>=vss_set))
                            if STABLE(g)==1 
                                for k=1:size(S_rov_aut_all,3)
                                    S_rov_aut_all2(:,:,k,g)=[kron(ones(size(incomevals)),S_rov_aut_all(:,:,k,g))];
                                    idio_dist_aut_all2(:,:,k,g)=[kron(ones(size(incomevals)),idio_dist_aut_all(:,:,k,g))];
                                end
                            end
                        end
                        S_rov_aut_all=S_rov_aut_all2;
                        idio_dist_aut_all=idio_dist_aut_all2;
                        clear S_rov_aut_all2;
                        clear idio_dist_aut_all2;
                    end
                    
                    P=kron(Q1,P_rov);
                    
                    Y=S(:,1)+ScaleRov*S(:,2);			% possible aggregate income outcome
                    state=length(S);                    % number of states
                    agent=2;                            % number of 'agents' in the HH vs. RoV approximation equals 2
                    STATE(vss)=state;
                    
                    
                    % ===================================
                    % Compute set of autarky values
                    % ===================================
                    
                    for k=1:agent
                        caut(:,k)=S(:,k);
                    end
                    waut=zeros(state,agent);
                    waut_avg_lag=zeros(state,agent);
                    waut_sb=zeros(state,agent);
                    % this is just the PDV of utility from consumption of
                    % income caut
                    eeewaut(:,1)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,caut(:,1));
                    eeewaut(:,2)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,caut(:,2));
                    
                    %%this is needed to compute the deviation values if we
                    %%want to take into account the previous history and
                    %%not always deviate with equal sharing
                     gammagrid=[lambda vss/ScaleRov*(1-lambda)];
                        gammagrid=[gammagrid ones(m,1) gammagrid(:,2)./gammagrid(:,1)];
                        gammagridold=gammagrid;
                        clear gammagrid
                        gammagrid=[linspace(1.2*max(S(:,1)),min(S(:,1)),m)',linspace(min(S(:,2)),1.2*max(S(:,2)),m)',ones(m,1),linspace(1.2*max(S(:,1)),min(S(:,1)),m)'./linspace(min(S(:,2)),1.2*max(S(:,2)),m)'];
                        
                   
                        % calculate the autarky value for ROV from average utility
                        % of its members in village of size vss
                        if vss>1 && Rov_sb==1
                            for j=1:state/3
                                % loop through individuals in the ROV
                                if lagind==1
                                    for jj=1:length(idio_dist(j,:))
                                        if idio_dist(j,jj)>0
                                            % distribution of incomes in ROV without one
                                            % individual
                                            idio_dist_lag=idio_dist(j,:);
                                            idio_dist_lag(1,jj)=idio_dist_lag(1,jj)-1;
                                            % make a temporary new 'state' of individual
                                            % and rov
                                            S_lag_temp(jj,2)=idio_dist_lag*incomevals;
                                            S_lag_temp(jj,1)=incomevals(jj);
                                            % find the corresponding state in a village of
                                            % size vss-1
                                            S_lag_2_temp=( ones(length(S_lag)/3,1)*idio_dist_lag==IDIO_DIST_ROV_KEEP(1:length(S_lag)/3,:,vss-1)) ;
                                            S_lag_2=sum(S_lag_2_temp,2)==3;
                                            S_lag_2=[S_lag_2;S_lag_2;S_lag_2];
                                            S_lag_1=(S_lag_temp(jj,1)==S_lag(:,1));
                                            S_lag_ind=find(S_lag_1.*S_lag_2==1);
                                            % construct utility of individual jj
                                            if one_per_aut==1
                                                waut_temp(jj,1)=utility(rho,S_lag_temp(jj,1)*(1-Pun_fac))+beta*P_lag(S_lag_ind,:)*waut_lag(:,1);
                                            else
                                                waut_temp(jj,1)=waut_lag(S_lag_ind,1);
                                            end
                                        else
                                            waut_temp(jj,1)=0;
                                        end
                                    end
                                    
                                    % utility of individuals is average of its members
                                    waut_rov(j,1)=idio_dist(j,:)*waut_temp/sum(idio_dist(j,:));
                                    
                                   
                                elseif lagind>1
                                    waut_rov(j,1)=100000;
                                end
                            end
                            waut_rov_sb=[waut_rov;waut_rov;waut_rov];
                        end

                        % in first iteration, or for individual deviations, waut is just consumption of income
                        if vss==1 || coaldev==0
                            for j=1:state
                                if one_per_aut==1
                                waut(j,1)=utility(rho,caut(j,1)*(1-Pun_fac))+beta*P(j,:)*eeewaut(:,1);
                                % autarky values for the rest of village are always the same - perfect
                                % risk-sharing among the rest of village
                                if Rov_sb==0 || vss==1
                                    waut(j,2)=utility(rho,caut(j,2)*(1-Pun_fac))+beta*P(j,:)*eeewaut(:,2);
                                else
                                    waut(j,2)=utility(rho,caut(j,2)*(1-Pun_fac))+beta*P(j,:)*waut_rov_sb;
                                end   
				else
					waut(j,1)=eeewaut(j,1);
                                		% autarky values for the rest of village are always the same - perfect
                               		 % risk-sharing among the rest of village
                                		if Rov_sb==0 || vss==1
                                	    waut(j,2)=eeewaut(j,2);
                                	else
                                    waut(j,2)=waut_rov_sb(j,1);
                                	end   

				end
                            end
                        else
                            % this calculates autarky values for the coalitional
                            % deviations case
                            
                            
                            for j=1:state
                                % indicator for the 'autarky' state in the state
                                % space of the previous iteration S_lag
                                S_lag_2=(abs(S_rov_aut(j)-S_lag(:,2))<1e-8);%(S_rov_aut(j)==S_lag(:,2)) ;
                                S_lag_1=(abs(S(j,1)-S_lag(:,1))<1e-8);%(S(j,1)==S_lag(:,1));
                                S_lag_ind=find(S_lag_1.*S_lag_2==1);
                                % for the individual, autarky value is utility of
                                % current consumption plus relevant trans
                                % probability times waut_lag, the vector of values
                                % from the next smaller insurance group
                                if one_per_aut==1
                                    waut(j,1)=utility(rho,caut(j,1)*(1-Pun_fac))+beta*P_lag(S_lag_ind,:)*waut_lag(:,1);
                                    % for rest of village, the autarky value is just
                                    % expected disc utility of current consumption
                                    if Rov_sb==0
                                        waut(j,2)=utility(rho,caut(j,2)*(1-Pun_fac))+beta*P(j,:)*eeewaut(:,2);
                                    else
                                        waut(j,2)=utility(rho,caut(j,2)*(1-Pun_fac))+beta*P(j,:)*waut_rov_sb;
                                    end
                                else
                                    waut(j,1)=waut_lag(S_lag_ind,1);
                                    if Rov_sb==0
                                        waut(j,2)=eeewaut(j,2);
                                    else
                                        waut(j,2)=waut_rov_sb(j,1);
                                    end
                                end
                                
                                %%%%%Here comes the big code that cycles through all
                                %%%%%potential coalitions of people
                                for g=1:size(S_rov_aut_all,4)
                                    if STABLE(g)==1
                                        for k=1:size(S_rov_aut_all,3)
                                            S_lag_2=(S_rov_aut_all(j,1,k,g)==S_LAG(:,2,g)) ;
                                            S_lag_1=(S(j,1)==S_LAG(:,1,g));
                                            S_lag_ind=find(S_lag_1.*S_lag_2==1);
                                            
                                            agentstate=[];
                                            for kk=1:length(incomevals)
                                                if idio_dist_aut_all(j,kk,k,g)>0
                                                    agentstate=[agentstate kk];
                                                end
                                            end
                                            
                                            if size(S_lag_ind)>0
                                                for kk=1:length(agentstate) %%%note, you only need one entry per income state, not one entry per agent
                                                    
                                                    if one_per_aut==1
                                                        waut_check_ind(j,1,k,g)=utility(rho,S(j,1))+beta*P_LAG(S_lag_ind,find((P_LAG(S_lag_ind,:,g)>0)),g)*WAUT_LAG(find((P_LAG(S_lag_ind,:,g)>0)),1,g);
                                                        waut_check(j,kk,k,g)=utility(rho,incomevals(agentstate(kk)))+beta*P_LAG(S_lag_ind,find((P_LAG(S_lag_ind,:,g)>0)),g)*WAUT_LAG(find((P_LAG(S_lag_ind,:,g)>0)),2,g);
                                                    else
                                                        waut_check_ind(j,1,k,g)=WAUT_LAG(S_lag_ind,1,g);
                                                        waut_check(j,kk,k,g)=WAUT_LAG(S_lag_ind,2,g);
                                                    
                                                        waut_check_ind_lambda(:,j,1,k,g)=WAUT_LAG_LAMBDA(:,S_lag_ind,1,g);
                                                        waut_check_lambda(:,j,kk,k,g)=WAUT_LAG_LAMBDA(:,S_lag_ind,2,g);
                                                    
                                                    end
                                                    S_check(j,kk,k,g)=agentstate(kk);  
                                                    
                                                end
                                            end
                                            
                                        end
                                    end
                                end
                            end
                            
                            
                        end
                        
                        % ====================================================
                        % Compute constrained insurance
                        % ====================================================
                        
                        %%Note that if you want to remember the previous
                        %%Pareto weight in the deviations, then you already
                        %%need to define the gammagrid earlier above. 
                       
                        % ====== make grids of weights, consumption, and Lagrange multipliers for updating ======
                        n=length(gammagrid);
                        w=zeros(n,state,agent);
                        gammar=zeros(n,state,agent);
                        c=zeros(n,state,agent);
                        phi=zeros(n,state,agent);
                        
                        % 'first-best' consumption given aggregate resources and planner weights
                        for j=1:state
                            c(:,j,1)=Y(j)./(1+ScaleRov*gammagrid(:,4));
                            c(:,j,2)=(Y(j)-c(:,j,1))/ScaleRov;
                        end
                        cfirstbest=c;
                        %  relative planner weights
                        for j=1:state
                            for k=1:agent
                                gammar(:,j,k)=gammagrid(:,agent+k);
                            end
                        end
                        % Compute the first best value functions without binding
                        % constraints
                        w=zeros(n,state,agent); dww=1;
                        wold=1/(1-beta)*utility(rho,cfirstbest(:,:,:));
                        while dww>0.00001
                            w(:,:,1)=utility(rho,cfirstbest(:,:,1))+beta*wold(:,:,1)*P';
                            w(:,:,2)=utility(rho,cfirstbest(:,:,2))+beta*wold(:,:,2)*P';
                            dww=max(max(max(abs(w-wold))));
                            wold=w;
                        end
                        wfirstbest=w;
                        % continuation utility as function of planner weight and state of the
                        % economy
                        ewfirstbest(:,:,1)=wfirstbest(:,:,1)*P';
                        ewfirstbest(:,:,2)=wfirstbest(:,:,2)*P';
                        
                        % ====================================================
                        % Start the loop with the enforcement constraints
                        % ====================================================n
                        ewsb=ewfirstbest;
                        wsb=wfirstbest;
                        options=optimset('Display','off');
                        eps=0.0000001;t=0;itgap=0; GAP=2; gap=2; gaphist=gap; gapincrease=0;
                        csolold=cfirstbest;
                        if beta<0.85
                            gapcrit=gapcritset;
                            maxitgap=maxitgapset;
                        else
                            gapcrit=gapcritset*1
                            maxitgap=maxitgapset;
                        end
                        
                        %%idiosyncratic distribution in rov for all states 
                        while GAP >=gapcrit && itgap<=maxitgap
                            rr=[];
                            rr1=[];
                            rr2=[];
                            rrnone=[];
                            rrboth=[];
                            nonconv=zeros(state,3)*nan;
                            itgap=itgap+1;
                            coalnonsust=nan(state,2);
                            count=0; countgrid=[ones(n,state)];
                            for j=1:state
                                
                                
                                index=[1 j];
                                
                                % 1. constraint binding for agent 1
                                % (individual)
                                rr1=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)<waut(j,1) & utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2)) ;
                                if vss>1 & coaldev==1
                                for i=1:length(rr1)
                                    if equal==0
                                    devindex=rr1(i);
                                    elseif equal==1
                                    devindex=(n+1)/2;    
                                    end
                                    index=[i j devindex];
                                    rr1cd=[];
                                    for g=size(waut_check,4):-1:0 %%%check from the largest group down, 0 is an individual deviation
                                        
                                        if  STABLECHECK(g+1)==1 & countgrid(i,j)>0  %%%you need to add one here, because you've got the individual in this vector as well
                                            %%%check the largest possible coalition first:
                                            if g>0
                                                for k=1:size(waut_check,3)   
                                                    if waut_check(j,1,k,g)~=0
                                                        coalition=[utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1)-waut_check_ind_lambda(devindex,j,1,k,g)];
                                                        for kk=1:size(waut_check,2)
                                                            if waut_check(j,kk,k,g)~=0
                                                                coalition=[coalition utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2)-waut_check_lambda(devindex,j,kk,k,g)];
                                                            end
                                                        end
                                                        coalition;
                                                        
                                                        %%%Step 1: find the grid points for which you are now trying to find a solution
                                                        
                                                        %%%%do the others want to come along?
                                                        rrcheck=coalition<0;
                                                        rrcheck=prod(rrcheck,2);
                                                        rrcheck=find(rrcheck>0);
                                                        
                                                        if size(rrcheck>0) & countgrid(i,j)>0
                                                            %%%%you need to do the solutions here
                                                            bind=1;
                                                            
                                                            coalitionsolve=[utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)-waut_check_ind_lambda(devindex,j,1,k,g)];
                                                            for kk=1:size(waut_check,2)
                                                                if waut_check(j,kk,k,g)~=0
                                                                    coalitionsolve=[coalitionsolve utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)-waut_check_lambda(devindex,j,kk,k,g)];
                                                                end
                                                            end
                                                            rrchecksolve=coalitionsolve<0;
                                                            rrchecksolve=prod(rrchecksolve,2);
                                                            rrchecksolve=find(rrchecksolve>0);
                                                            
                                                            rr=rrchecksolve(end)+1; 
                                                            if rr>length(gammagrid)
                                                                coalnonsust(j,1)=-2;
                                                                output=[nan,nan,-50];
                                                            else
                                                                
                                                                clear x0
                                                                
                                                                x0(1) = cfirstbest(i,j,1);  %%redundant because we no longer need an initial guess
                                                                %T vectors of consumption and utilty that meet
                                                                %participation constraints
                                                                [x, output]=agent1maxrho_rov_lambda_constant(x0,index,ewsb, waut_check_ind_lambda(:,:,1,k,g), waut(:,:) ,bind,Y,n,beta,rho,ScaleRov,cfirstbest,rr);
                                                                
                                                                %T 29_Mar This keeps track of nonconvergence, but
                                                                %allows the programme to continue - sometimes it
                                                                %converges in a subsequent iteration
                                                            end
                                                            if output(3) <=0
                                                                nonconv(j,1)=output(3);
                                                                csol(i,j,1)=caut(j,1);
                                                                csol(i,j,2)=caut(j,2);
                                                                gammar(i,j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                                                wsbnewsol(i,j,1)=waut(j,1);
                                                                wsbnewsol(i,j,2)=waut(j,2);
                                                                phisol(i,j,1)=gammagrid(i,4).^rho./gammar(i,j,2).^rho-1;
                                                                phisol(i,j,2:agent)=0	;
                                                            else
                                                                csol(i,j,1)=x(1);
                                                                gammar(i,j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                                                csol(i,j,2)=(Y(j)-x(1))/ScaleRov;
                                                                wsbnewsol(i,j,1)=output(1);
                                                                wsbnewsol(i,j,2)=output(2);
                                                                phisol(i,j,1)=gammagrid(i,4).^rho./gammar(i,j,2).^rho-1;
                                                                phisol(i,j,2:agent)=0	;
                                                            end
                                                            countgrid(i,j)=0;  %%%these grid points have been dealt with
                                                
                                                        end
                                                    end
                                                end
                                                
                                                  
                                            elseif g==0 & countgrid(i,j)>0%%%only an individual deviation
                                                coalition=[utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1)-eeewaut(j,1)];
                                                
                                                %%%%do the others want to come along?
                                                rrcheck=coalition<0;
                                                rrcheck=prod(rrcheck,2);
                                                rrcheck=find(rrcheck>0);
                                                if size(rrcheck>0)
                                                    %%%%you need to do the solutions here
                                                   bind=1;
                                                    coalitionsolve=[utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)-eeewaut(j,1)];
                                                    %%%%do the others want to come along?
                                                    rrchecksolve=coalitionsolve<0;
                                                    rrchecksolve=prod(rrchecksolve,2);
                                                    rrchecksolve=find(rrchecksolve>0);
                                                    
                                                   
                                                    rr=rrchecksolve(end)+1;
                                                    
                                                    if rr>length(gammagrid)
                                                        coalnonsust(j,1)=-2;
                                                        output=[nan,nan,-50];
                                                    else
                                                        clear x0
                                                        x0(1) = cfirstbest(rr,j,1);  %%redundant because we no longer need an initial guess
                                                        %T vectors of consumption and utilty that meet
                                                        %participation constraints
                                                        [x, output]=agent1maxrho_rov_06_Mar(x0,index,ewsb,[eeewaut],bind,Y,n,beta,rho,ScaleRov,cfirstbest,rr,coaldev);
                                                        
                                                        %T 29_Mar This keeps track of nonconvergence, but
                                                        %allows the programme to continue - sometimes it
                                                        %converges in a subsequent iteration
                                                    end
                                                    if output(3) <=0
                                                        nonconv(j,1)=output(3);
                                                        csol(i,j,1)=caut(j,1);
                                                        csol(i,j,2)=caut(j,2);
                                                        gammar(i,j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                                        wsbnewsol(i,j,1)=waut(j,1);
                                                        wsbnewsol(i,j,2)=waut(j,2);
                                                        phisol(i,j,1)=gammagrid(i,4).^rho./gammar(i,j,2).^rho-1;
                                                        phisol(i,j,2:agent)=0	;
                                                    else
                                                        csol(i,j,1)=x(1);
                                                        gammar(i,j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                                        csol(i,j,2)=(Y(j)-x(1))/ScaleRov;
                                                        wsbnewsol(i,j,1)=output(1);
                                                        wsbnewsol(i,j,2)=output(2);
                                                        phisol(i,j,1)=gammagrid(i,4).^rho./gammar(i,j,2).^rho-1;
                                                        phisol(i,j,2:agent)=0	;
                                                    end
                                                countgrid(i,j)=-1;  %%%because these grid points have been dealt with
                                                else %%%if there is no binding constraint at this grid point 
                                                csol(i,j,1:agent)=cfirstbest(i,j,1:agent);
                                                gammar(i,j,1:agent)=gammagrid(i,agent+1:2*agent);
                                                phisol(i,j,1:agent)=0;
                                    
                                                wsbnewsol(i,j,1)=utility(rho,cfirstbest(i,j,1))+beta*ewsb(i,j,1);
                                                wsbnewsol(i,j,2)=utility(rho,cfirstbest(i,j,2))+beta*ewsb(i,j,2);
                                                countgrid(i,j)=-2; 
                                                end
                                               
                                            end
                                        end
                                    end
                                end
                                %and this is the remainder;
                                    RRCHECK=countgrid(rr1,j)<=0;
                                    rr1(RRCHECK)=[];
                                    rrnoneadd=rr1;
                                else
                                    %%%%in first group size iteration, i.e. for vss=1
                                    rr1cd=rr1;
                                    rr1id=[];
                                    rrnoneadd=[];
                                end

                                %%%%leave this here for vss=1, where there are no
                                %%%%possible joint deviations
                                if size(rr1cd)>0
                                    bind=1;
                                    
                                    rr=rr1cd(end)+1;
                                    
                                    if isempty(rr)==1
                                        coalnonsust(j,1)=-2;
                                        output=[nan,nan,-50];
                                    else
                                        if size(rr1cd)>0
                                            clear x0
                                            x0(1) = cfirstbest(rr1cd(1),j,1);  %%redundant because we no longer need an initial guess
                                            %T vectors of consumption and utilty that meet
                                            %participation constraints
                                            [x, output]=agent1maxrho_rov_06_Mar(x0,index,ewsb,waut,bind,Y,n,beta,rho,ScaleRov,cfirstbest,rr,coaldev);
                                            
                                            %T 29_Mar This keeps track of nonconvergence, but
                                            %allows the programme to continue - sometimes it
                                            %converges in a subsequent iteration
                                        end
                                        
                                    end
                                    if output(3) <=0
                                        nonconv(j,1)=output(3);
                                        csol(rr1cd,j,1)=ones(size(rr1cd,1),1)*caut(j,1);
                                        csol(rr1cd,j,2)=ones(size(rr1cd,1),1)*caut(j,2);
                                        gammar(rr1cd,j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                        wsbnewsol(rr1cd,j,1)=ones(size(rr1cd,1),1)*waut(j,1);
                                        wsbnewsol(rr1cd,j,2)=ones(size(rr1cd,1),1)*waut(j,2);
                                        phisol(rr1cd,j,1)=gammagrid(rr1cd,4).^rho./gammar(rr1cd,j,2).^rho-1;
                                        phisol(rr1cd,j,2:agent)=0	;
                                    else
                                        csol(rr1cd,j,1)=x(1);
                                        gammar(rr1cd,j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                        csol(rr1cd,j,2)=(Y(j)-x(1))/ScaleRov;
                                        wsbnewsol(rr1cd,j,1)=output(1);
                                        wsbnewsol(rr1cd,j,2)=output(2);
                                        phisol(rr1cd,j,1)=gammagrid(rr1cd,4).^rho./gammar(rr1cd,j,2).^rho-1;
                                        phisol(rr1cd,j,2:agent)=0	;
                                    end
                                end
                                %this becomes relevant for vss>1 and mops
                                %up grid points that were not dealt with
                                %above.
                                if size(rrnoneadd)>0
                                    
                                    csol(rrnoneadd,j,1:agent)=cfirstbest(rrnoneadd,j,1:agent);
                                    gammar(rrnoneadd,j,1:agent)=gammagrid(rrnoneadd,agent+1:2*agent);
                                    phisol(rrnoneadd,j,1:agent)=0;
                                    
                                    wsbnewsol(rrnoneadd,j,1)=utility(rho,cfirstbest(rrnoneadd,j,1))+beta*ewsb(rrnoneadd,j,1);
                                    wsbnewsol(rrnoneadd,j,2)=utility(rho,cfirstbest(rrnoneadd,j,2))+beta*ewsb(rrnoneadd,j,2);
                                end
                                
                                % 2. Constraint binding for agent 2 (RoV)
                                rr2=find(utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)<waut(j,2)) ;
                                
                                if size(rr2)>0
                                    bind=2; %%i.e. agent 2 has the binding constraint
                                    rr=find(utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2));
                                    if isempty(rr)==1
                                        coalnonsust(j,2)=-2;
                                        output=[nan,nan,-50];
                                    else
                                        clear x0
                                        x0(1) = cfirstbest(rr(end),j,1); %%redundant because we no longer need an initial guess
                                        %%%%note that I'm now checking whether
                                        %%%%individual is happy with this allocation
                                        %%%%relative to autarky, given that rov gets
                                        %%%%waut(j,2) inside contract, we know that
                                        %%%%there is nothing that the
                                        [x, output]=agent1maxrho_rov_06_Mar(x0,index,ewsb,[eeewaut(:,1) waut(:,2)],bind,Y,n,beta,rho,ScaleRov,cfirstbest,rr,coaldev);
                                        %%%%check that agent 1 does not have a binding
                                        %%%%constraint at this new solution, and note that
                                        %%%%the rov must never have less than waut(j,2),
                                        %%%%so if individual has a binding constraint that
                                        %%%%means rov would have to decrease its
                                        %%%%consumption, not possible, so you do not need
                                        %%%%to solve, just check!
                                    end
                                    % output(3) equals 1 if constraints
                                    % are satisfied, but -2 if not
                                    if output(3)<=0
                                        nonconv(j,2)=output(3)*100;
                                        csol(rr2,j,1)=ones(size(size(rr2,1),1),1)*caut(j,1);
                                        gammar(rr2,j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                        csol(rr2,j,2)=ones(size(size(rr2,1),1),1)*caut(j,2);
                                        wsbnewsol(rr2,j,1)=ones(size(rr2))*waut(j,1);
                                        wsbnewsol(rr2,j,2)=ones(size(rr2))*waut(j,2);
                                        phisol(rr2,j,2)=gammar(rr2,j,2).^rho./gammagrid(rr2,4).^rho-1;
                                        phisol(rr2,j,1)=0;
                                    else
                                        %%%%first check that agent 1 (in conjunction with the rest of the village)
                                        %%%%truly does not have binding constraints
                                    if vss>1  & coaldev==1 % only need to check this once coalitional deviations become relevant
                                        for i=1:length(rr2)
                                        if equal==0
                                        devindex=rr2(i);%i;%(n+1)/2;%i;%(n+1)/2;
                                        elseif equal==1
                                        devindex=(n+1)/2;    
                                        end
                                            for g=size(waut_check,4):-1:1 %%%check from the largest group down, 0 is an individual deviation
                                                 
                                                 if  STABLECHECK(g+1)==1   %%%you need to add one here, because you've got the individual in this vector as well
                                                     %%%check the largest possible coalition first:
                                                     for k=1:size(waut_check,3)   %%%
                                                         if waut_check(j,1,k,g)~=0
                                                             coalition=[output(1)-waut_check_ind_lambda(devindex,j,1,k,g)];
                                                             for kk=1:size(waut_check,2)
                                                                 if waut_check(j,kk,k,g)~=0
                                                                     coalition=[coalition output(2)-waut_check_lambda(devindex,j,kk,k,g)];
                                                                 end
                                                             end
                                                             coalition;
                                                             %%%Step 1: find the grid points for which you are now trying to find a solution
                                                             %%%%do the others want to come along?
                                                             rrcheck=coalition<0;
                                                             rrcheck=prod(rrcheck,2);
                                                             rrcheck=find(rrcheck>0);
                                                             if size(rrcheck)>0
                                                                 i
                                                                 output(3)=-2;  %%%yes, so if any of these have a binding constraint, then output(3) is set equal to -2;
                                                             end
                                                         end
                                                     end
                                                 end
                                             end
                                         end
                                    end
                                        if output(3)>0
                                            csol(rr2,j,1)=x(1);
                                            gammar(rr2,j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                            csol(rr2,j,2)=(Y(j)-x(1))/ScaleRov;
                                            wsbnewsol(rr2,j,1)=output(1);
                                            wsbnewsol(rr2,j,2)=output(2);
                                            phisol(rr2,j,2)=gammar(rr2,j,2).^rho./gammagrid(rr2,4).^rho-1;
                                            phisol(rr2,j,1)=0;
                                        else
                                            nonconv(j,2)=output(3)*100;
                                            csol(rr2,j,1)=ones(size(size(rr2,1),1),1)*caut(j,1);
                                            gammar(rr2,j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                            csol(rr2,j,2)=ones(size(size(rr2,1),1),1)*caut(j,2);
                                            wsbnewsol(rr2,j,1)=ones(size(rr2))*waut(j,1);
                                            wsbnewsol(rr2,j,2)=ones(size(rr2))*waut(j,2);
                                            phisol(rr2,j,2)=gammar(rr2,j,2).^rho./gammagrid(rr2,4).^rho-1;
                                            phisol(rr2,j,1)=0;
                                        end
                                    end
                                end
                                
                                % 3. Constraint all slack
                                rrnone=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)>=waut(j,1) & utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2)) ;
                                if size(rrnone)>0
                                    csol(rrnone,j,1:agent)=cfirstbest(rrnone,j,1:agent);
                                    gammar(rrnone,j,1:agent)=gammagrid(rrnone,agent+1:2*agent);
                                    phisol(rrnone,j,1:agent)=0;
                                    wsbnewsol(rrnone,j,1)=utility(rho,cfirstbest(rrnone,j,1))+beta*ewsb(rrnone,j,1);
                                    wsbnewsol(rrnone,j,2)=utility(rho,cfirstbest(rrnone,j,2))+beta*ewsb(rrnone,j,2);
                                end
                            end
                            gapold=max(norm(wsbnewsol(:,:,1)-wsb(:,:,1)),norm(wsbnewsol(:,:,2)-wsb(:,:,2)));
                            gap1=max(max(max(abs(invutility(beta,rho,wsbnewsol(:,:,1))-invutility(beta,rho,wsb(:,:,1)))./invutility(beta,rho,wsb(:,:,1)))),max(max(abs(invutility(beta,rho,wsbnewsol(:,:,2))-invutility(beta,rho,wsb(:,:,2)))./invutility(beta,rho,wsb(:,:,2)))))/(1-beta);
                            gap2=max(max(max(abs(csol(:,:,1)-csolold(:,:,1))./csol(:,:,1))),max(max(abs(csol(:,:,2)-csolold(:,:,2))./csol(:,:,2))));
                            gap=max(gap1,gap2);
                            gaphist(itgap+1)=gap;
                            GAP=max(gaphist(max(1,length(gaphist)-convlength):end))
                            disp('gap,gapold')
                            disp([gap,gapold])
                            nonconvind(Qbeta,Qrho,vss,incproc)=mean(nonconv(find(~isnan(nonconv))));
                            Coalnonsust(Qbeta,Qrho,vss,incproc)=mean(coalnonsust(find(~isnan(coalnonsust))));
                            
                           
                            if  (min(min(nonconv))<=0 ||min(min(coalnonsust))<0) & itgap>=3 & coaldev==1
                                itgap
                                disp('Autarky or coalition not sustainable')
                                [nonconv,S]
                                for j=1:state
                                    csol(:,j,1)=caut(j,1);
                                    csol(:,j,2)=caut(j,2);
                                end
                                break
                            end
                            csolold=csol;
                            wsb=wsbnewsol;
                            ewsb(:,:,1)=(P*wsb(:,:,1)')';
                            ewsb(:,:,2)=(P*wsb(:,:,2)')';
                        end
                        CRITGAP(vss)=gaphist(itgap-1);
                        if itgap>=maxitgap
                            nonconvind(Qbeta,Qrho,vss,incproc)=gap;
                        end
                        if gapincrease==1
                            nonconvind(Qbeta,Qrho,vss,incproc)=-gap;
                        end
                        
                        % =============================================================
                        % calculate values of outside option for coalitional deviations
                        % =============================================================
                        
                        if coaldev==1 && (nonconvind(Qbeta,Qrho,vss,incproc)<=0 || Coalnonsust(Qbeta,Qrho,vss,incproc)<=0) && itgap>=3 && vss>1
                            % if vss is unsustainable, the outside option for vss+1 remains that
                            % for vss
                            waut_lag(:,1)=waut_lag(:,1);
                            waut_lag(:,2)=waut_lag(:,2);
                            P_lag=P_lag;
                            S_lag=S_lag;
                            % if the current insurance group is not sustainable, it
                            % cannot be the outside option for the next bigger group,
                            % so increase the 'lagind' indicator by one
                            lagind=lagind+1;
                            gap=-1;
                            STABLE(vss)=0;
                            STABLECHECK=[1 STABLE];
                            % ... or if there is no previous, use just conditional exp of
                            % autarky consumption as calculated above
                        elseif (nonconvind(Qbeta,Qrho,vss,incproc)<=0 || Coalnonsust(Qbeta,Qrho,vss,incproc)<=0)  && itgap>=3 && (vss==1 || coaldev==0)
                            waut_lag(:,1)=eeewaut(:,1);
                            waut_lag(:,2)=eeewaut(:,2);
                            P_lag=P;
                            S_lag=S;
                            lagind=1;
                            gap=-1;
                            
                            S_rov_KEEP(1:state/3,vss)=S_rov;
                            S_KEEP(1:state,:,vss)=S;
                            PHI_KEEP(:,1:state,:,vss)=phisol;
                            STABLE(vss)=1;
                            STABLECHECK=[1 STABLE];
                            IDIO_DIST_ROV_KEEP(1:state/3,1:3,vss)=idio_dist;
                        else 
                            lagind=1;
                            % when ins group of size vss is sustainable, its value
                            % becomes outside option for next bigger group size
                            for j=1:state
                                % outside option for next larger insurance group is
                                % value with equal weights vss ...
                                eewaut(j,1)=wsb((n+1)/2,j,1);
                                eewaut(j,2)=wsb((n+1)/2,j,2);
                                for i=1:n
                                eewaut_lambda(i,j,1)=wsb(i,j,1);
                                eewaut_lambda(i,j,2)=wsb(i,j,2);
                                end
                            end
                            waut_lag=eewaut;
                            waut_lag_lambda=eewaut_lambda;
                            S_lag=S;
                            P_lag=P;
                            S_rov_KEEP(1:state/3,vss)=S_rov;
                            S_KEEP(1:state,:,vss)=S;
                            PHI_KEEP(:,1:state,:,vss)=phisol;
                            WAUT_LAG(1:state,1:2,vss)=waut_lag;
                            WAUT_LAG_LAMBDA(1:n,1:state,1:2,vss)=waut_lag_lambda;
                            P_LAG(1:state,1:state,vss)=P_lag;
                            S_LAG(1:state,:,vss)=S_lag;
                            STABLE(vss)=1;
                            STABLECHECK=[1 STABLE]; %%%this includes stability for the individual as well, makes the loop above easier
                            IDIO_DIST_ROV_KEEP(1:state/3,1:3,vss)=idio_dist;
                        end
                        if itgap>=maxitgap
                            nonconvind(Qbeta,Qrho,vss,incproc)=gap;
                        end
                        
                    
                
                
                
            end
            
           
            %==========================================================
            % Simulate the economy
            % =========================================================
            % This lists all possible combinations of individual income
            % values (IDIO_DIST), and the aggregate states (pair of individual income and ROV-income) associated
            % with every individual income value

            if coaldev==1                       
                stable=(find(STABLE==1));
                vss_set_sim=stable;
            elseif coaldev==0                             
                vss_set_sim=vss_set;
            elseif coaldev==2
                vss_set_sim=vss_set;
            elseif coaldev==3
                vss_set_sim=stable';
            end

            for vss=[vss_set_sim]
                if coaldev<2
                    state=STATE(vss);
                    if rov_average==1
                        ScaleRov=vss;
                    else
                        ScaleRov=1;
                    end
                    S=S_KEEP(1:state,:,vss);
                    S_rov=S_rov_KEEP(1:state/3,:,vss);
                    idio_dist=IDIO_DIST_ROV_KEEP(1:state/3,1:3,vss);
                    phisol=PHI_KEEP(:,1:state,:,vss);
                    IDIO_DIST=[ones(state/3,1)*[1 0 0] + idio_dist; ones(state/3,1)*[0 1 0] + idio_dist; ones(state/3,1)*[0 0 1]+idio_dist];
                    distribution=size(IDIO_DIST,1);
                end
                
                Nvill=ceil((vss_high_inddev+1)/(vss+1))
                var_dyagglog=zeros(Nvill*(vss+1),vss_high);var_dyvillagglog=zeros(Nvill*(vss+1),vss_high);vardaggylog=zeros(Nvill*(vss+1),vss_high);betayc=zeros(2,Nvill*(vss+1));vardcshare=zeros(Nvill*(vss+1),vss_high);vardcshare_dylow=zeros(Nvill*(vss+1),vss_high);vardcshare_dyhigh=zeros(Nvill*(vss+1),vss_high);varclog=zeros(Nvill*(vss+1),vss_high);vardclog=zeros(Nvill*(vss+1),vss_high);varylog=zeros(Nvill*(vss+1),vss_high);vardylog=zeros(Nvill*(vss+1),vss_high);varclog_cond=zeros(Nvill*(vss+1),vss_high);varylog_cond=zeros(Nvill*(vss+1),vss_high);
                varc_yhigh_cond=zeros(Nvill*(vss+1),vss_high);varc_ylow_cond=zeros(Nvill*(vss+1),vss_high);vardclog_cond=zeros(Nvill*(vss+1),vss_high);vardylog_cond=zeros(Nvill*(vss+1),vss_high); varc_yhigh=zeros(Nvill*(vss+1),vss_high);
                vardc_dypos_cond=zeros(Nvill*(vss+1),vss_high);
                vardc_dyneg_cond=zeros(Nvill*(vss+1),vss_high); varc_yhigh=zeros(Nvill*(vss+1),vss_high); varc_ylow=zeros(Nvill*(vss+1),vss_high);vardc_dypos=zeros(Nvill*(vss+1),vss_high); vardc_dyneg=zeros(Nvill*(vss+1),vss_high);
                R2_dcdy=zeros(1,Nvill*(vss+1));betadcdy_agg=zeros(3,Nvill*(vss+1));   betadcdy=zeros(2,vss_high);  betadcdy_high=zeros(1,Nvill*(vss+1));betadcdy_low=zeros(1,Nvill*(vss+1));relcvar=zeros(Nvill*(vss+1),vss_high); reldcvar=zeros(Nvill*(vss+1),vss_high); reldc0var=zeros(Nvill*(vss+1),vss_high); relcvar_cond=zeros(Nvill*(vss+1),vss_high); reldcvar_cond=zeros(Nvill*(vss+1),vss_high);
                beta_cond=zeros(2,Nvill*(vss+1),vss_high); beta_ylow_cond=zeros(4,Nvill*(vss+1),vss_high); beta_yhigh_cond=zeros(4,Nvill*(vss+1),vss_high) ;
                beta_agg=zeros(2,Nvill*(vss+1),vss_high); beta_yhigh_agg=zeros(2,Nvill*(vss+1),vss_high); beta_ylow_agg=zeros(2,Nvill*(vss+1),vss_high);
                
                %=============================================
                % Start simulation
                % ============================================
                % number of insurance groups equals village size
                % vss_high+1 divided by size of ins group vss+1
                
                Npanel=Nvill*(vss+1);
                csim0=nan(Nvill,vss+1,T);
                income0=nan(Nvill,vss+1,T);
                income0ind=nan(Nvill,vss+1,T);
                phisim0=nan(Nvill,vss+1,T);
                Yvill=nan(1,T);
                Yagg=nan(Nvill,T);
                Ssim=nan(Nvill,vss+1,T);
                gammarfun=nan(Nvill,vss+1,T);
                phinew=nan(Nvill,vss+1,T);
                indSjj=nan(T,vss+1);
                disp('simulating') 
                disp('Reminder: now simulating for village,coaldev,incproc')
                disp([village coaldev incproc])
                disp('now simulating for parameters')
                disp([Qbeta,Qrho])
                disp('vss')
                disp(vss)
                tic
                %initial state
                t=1;
                % initial values of income and consumption (set =
                % income)
                rng(3,'twister');
                shockinittemp=randi(3,Nvill*(vss+1),1);
                shockinit=reshape(shockinittemp,Nvill,vss+1);
                income0ind(:,:,1)=shockinit;
                income0(:,:,1)=incomevals(income0ind(:,:,1));
                csim0(:,:,1)=income0(:,:,1);
                Yagg(:,t)=sum(income0(:,:,t),2);
                IDIO_DIST_sim=[];
                
                while t<T
                    t=t+1;
                    %%%%simulate next period income taking into account
                    %%%% persistence across time
                    rng(t,'twister');
                    shocktemp=rand(Nvill*(vss+1),1);
                    shock=reshape(shocktemp,Nvill,vss+1);
                    for k=1:vss+1
                        rr=(find(shock(:,k)<=cumQ(income0ind(:,k,t-1),1)));
                        income0ind(rr,k,t)=1;
                        rr=(find(cumQ(income0ind(:,k,t-1),1)<shock(:,k) & cumQ(income0ind(:,k,t-1),2)>=shock(:,k)));
                        income0ind(rr,k,t)=2;
                        rr=(find(cumQ(income0ind(:,k,t-1),2)<shock(:,k) & cumQ(income0ind(:,k,t-1),3)>=shock(:,k)));
                        income0ind(rr,k,t)=3;
                        income0(:,k,t)=incomevals(income0ind(:,k,t))  ;
                    end
                    Yagg(:,t)=sum(income0(:,:,t),2);
                end
                
                
                    t=1;
                    while t<T
                        t=t+1;
                        gammarfun(:,:,t)=(Yagg(:,t-1)*ones(1,vss+1)-csim0(:,:,t-1))./csim0(:,:,t-1)/ScaleRov;
                        gammarfun(:,:,t)=min(gammarfun(:,:,t),max(gammagrid(:,4))*ones(Nvill,vss+1));                        
                        for ivill=1:Nvill
                            IDIO_DIST_sim=find(((IDIO_DIST(:,1,1)==(sum(income0ind(ivill,:,t)==1))).*(IDIO_DIST(:,2,1)==(sum(income0ind(ivill,:,t)==2))).*(IDIO_DIST(:,3,1)==(sum(income0ind(ivill,:,t)==3))))==1)';
                            for k=1:vss+1
                                % choose for individual k the aggregate state
                                % that applies
                                indSjj(t,k)=IDIO_DIST_sim(S(IDIO_DIST_sim,1)==income0(ivill,k,t));
                                % check that aggregate and individual income0s are
                                % consistent with the states caculated above
                                if abs(Yagg(ivill,t)-(S(indSjj(t,k),1)+ScaleRov*S(indSjj(t,k),2)))>0.000001
                                    stop
                                end
                            end
                            Ssim(ivill,:,t)=indSjj(t,:);
                        end
                        for k=1:vss+1
                            BLA=interp1(gammagrid(:,4),phisol(:,Ssim(:,k,t),1),gammarfun(:,k,t));
                            phinew(:,k,t)=diag(BLA);
                        end
                        
                       
                        gammaraux=((1+phinew(:,:,t))./(1+phinew(:,1,t)*ones(1,vss+1))).^(1/rho).*csim0(:,:,t-1)./(csim0(:,1,t-1)*ones(1,vss+1));
                       
                        
                        csim0(:,:,t)=(gammaraux./(sum(gammaraux,2)*ones(1,vss+1))).*(Yagg(:,t)*ones(1,vss+1));
                       
                    end
                    AA=[];
                    BB=[];
                    CC=[];
                    for t=1:T
                        AA=[AA reshape(csim0(:,:,t)',Nvill*(vss+1),1)];
                        BB=[BB reshape(income0(:,:,t)',Nvill*(vss+1),1)];
                        CC=[CC reshape(income0ind(:,:,t)',Nvill*(vss+1),1)];
                    end
                    csim0=AA;
                    income0ind=CC;
                    income0=BB;
               
                disp('simulating takes')
                tsim(Qbeta,Qrho)=toc;            
                disp(tsim(Qbeta,Qrho))
                
                % =========================================================
                % calculate moments of consumption and income growth
                % =========================================================
                csim =[]; cshare_raw =[]; dcshare_raw =[]; cc =[]; ccc =[]; yy =[]; aggyy =[]; aggvillincome=[]; yyy_cond  =[]; dc_log_cond =[]; dy_log_cond =[]; ccc_cond=[];
                rng(3, 'twister');
                cerr=randn(size(csim0)); logcsim0_groupmeandev=[]; logcsim_groupmeandev=[]; logincome_groupmeandev=[];
                rng(4, 'twister');
                incerr=randn(size(csim0));   Momentsest1=nan(35,N_cerr); ind=[];
                % loop over size of cons measurement error
                for kkk=1:size(csim0,1)
                    tempvar(kkk)=var(log(csim0(kkk,:)));
                end
                Tempvar=mean(tempvar);
                for ic=1:N_cerr+1
                    if N_cerr>0
                        csim=exp(log(csim0)+((ic-1)/N_cerr*maxcerr/(1-(ic-1)/N_cerr*maxcerr)*Tempvar)^0.5*cerr);
                    else
                        csim=exp(log(csim0));
                    end
                    logcsim=log(csim);
                    logcsim0=log(csim0);
                    logcsim0_meandev=logcsim0-ones(size(csim0,1),1)*mean(logcsim0,1);
                    logcsim_meandev=logcsim-ones(size(csim0,1),1)*mean(logcsim,1);
		    
                    % add inc meas error back in
                    income=exp(log(income0)+varnu(incproc)^0.5*incerr);
                    logincome=log(income);
                    logincome0=log(income0);
                    logincome_meandev=logincome-ones(size(csim0,1),1)*mean(logincome,1);
                    income0_meandev=income0-ones(size(csim0,1),1)*mean(income0,1);
                    for ivill=1:Nvill
                        aggvillincome((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=ones(vss+1,1)*sum(income((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                        logaggvillincome=log(aggvillincome);
                        logcsim0_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logcsim0((ivill-1)*(vss+1)+1:ivill*(vss+1),:)-ones(vss+1,1)*mean(logcsim0((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                        logcsim_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logcsim((ivill-1)*(vss+1)+1:ivill*(vss+1),:)-ones(vss+1,1)*mean(logcsim((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                        logincome_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logincome((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)-ones(vss+1,1)*mean(logincome((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                    end
                    aggincome=sum(income,1);
                    % this chooses log or level consumption for
                    % moments
                    if momentclevinlog==1
                        ccc=logcsim(:,Tinit:T);
                        ccc_cond=logcsim_meandev(:,Tinit:T);
                        ccc_groupcond=logcsim_groupmeandev(:,Tinit:T);
                    else
                        ccc=(csim(:,Tinit:T));
                        ccc_cond=(csim_meandev(:,Tinit:T)); %NB this doesn't actually exist...
                    end
                    yyy=logincome(:,Tinit:T);
                    yyy_cond=logincome_meandev(:,Tinit:T);
                    yyy_groupcond=logincome_groupmeandev(:,Tinit:T);
                    cc=logcsim(:,Tinit:T)-logcsim(:,Tinit-1:T-1);
	             var_dc_t_fe(vss)=var(mean(cc,1));
                    cc_cond=logcsim_meandev(:,Tinit:T)-logcsim_meandev(:,Tinit-1:T-1);
                    cc_groupcond=logcsim_groupmeandev(:,Tinit:T)-logcsim_groupmeandev(:,Tinit-1:T-1);
                    yy=logincome(:,Tinit:T)-logincome(:,Tinit-1:T-1);
                    yy_cond=logincome_meandev(:,Tinit:T)-logincome_meandev(:,Tinit-1:T-1);
                    yy_groupcond=logincome_meandev(:,Tinit:T)-logincome_meandev(:,Tinit-1:T-1);
                    aggyy=log(sum(income(:,Tinit:T),1))-log(sum(income(:,Tinit-1:T-1),1));
                    aggvillyy=log(aggvillincome(:,Tinit:T))-log(aggvillincome(:,Tinit-1:T-1));
                    vardaggylog(1,vss)=var(log(aggincome(1,Tinit:T))-log(aggincome(1,Tinit-1:T-1)));
                    
                    % calculate moments for each individiual
                    for ind=1:size(csim,1)
                        dc_log_cond(ind,:)=(ccc_cond(ind,2:end))-(ccc_cond(ind,1:end-1));
                        dy_log_cond(ind,:)=(yyy_cond(ind,2:end))-(yyy_cond(ind,1:end-1));
                        dc_log_groupcond(ind,:)=(ccc_groupcond(ind,2:end))-(ccc_groupcond(ind,1:end-1));
                        dy_log_groupcond(ind,:)=(yyy_groupcond(ind,2:end))-(yyy_groupcond(ind,1:end-1));
                        
                        % conditioning on income increases, and income
                        % decreases
                        % NB: we allocate no change to decrease!!
                        rrhc =[]; rrlc =[]; rrl =[]; rrh=[];
                        rrhc=find(yy(ind,:)>0);
                        rrlc=find(yy(ind,:)<=0);
                        rrhc_cond=find(yy_cond(ind,:)>0);
                        rrlc_cond=find(yy_cond(ind,:)<=0);
                        rrhc_groupcond=find(yy_groupcond(ind,:)>0);
                        rrlc_groupcond=find(yy_groupcond(ind,:)<=0);
                        vardclog(ind,vss)=var(cc(ind,:));
                        vardylog(ind,vss)=var(yy(ind,:));
                        vardc_dypos(ind,vss)=var(cc(ind,rrhc));
                        vardc_dyneg(ind,vss)=var(cc(ind,rrlc));
                        
                        %%%conditional on village income
                        vardclog_cond(ind,vss)=var(dc_log_cond(ind,:));
                        vardylog_cond(ind,vss)=var(dy_log_cond(ind,:));
                        vardc_dypos_cond(ind,vss)=var(dc_log_cond(ind,dy_log_cond(ind,:)>0));
                        vardc_dyneg_cond(ind,vss)=var(dc_log_cond(ind,dy_log_cond(ind,:)<=0));
                        
                        %%%conditional on group income
                        vardclog_groupcond(ind,vss)=var(dc_log_groupcond(ind,:));
                        vardylog_groupcond(ind,vss)=var(dy_log_groupcond(ind,:));
                        vardc_dypos_groupcond(ind,vss)=var(dc_log_groupcond(ind,dy_log_groupcond(ind,:)>0));
                        vardc_dyneg_groupcond(ind,vss)=var(dc_log_groupcond(ind,dy_log_groupcond(ind,:)<=0));
                      
                        vardc_dypos(ind,vss)=var(cc(ind,rrhc));
                        vardc_dyneg(ind,vss)=var(cc(ind,rrlc));
                        
                        %%% Relative variances
                        reldcvar(ind,vss)=var(cc(ind,rrhc))-var(cc(ind,rrlc));
                        reldcvar_cond(ind,vss)=vardc_dypos_cond(ind,vss)-vardc_dyneg_cond(ind,vss);
                        reldcvar_groupcond(ind,vss)=vardc_dypos_groupcond(ind,vss)-vardc_dyneg_groupcond(ind,vss);
                        
                        %%% THE BETA COEFFICIENTS
                        %1) Unconditional
                        %risk-sharing whole sample
                        covtemp=cov(cc(ind,:),aggyy(1,:));
                        betadcdy_agg(ind,vss)=covtemp(2,1)/var(aggyy(1,:));
                        covtemp=cov(cc(ind,:),yy(ind,:));
                        betadcdy(ind,vss)=covtemp(2,1)/(var(yy(ind,:)));
                        [betadcdytemp,temp1,rtemp,temp2,statstemp]=regress(cc(ind,:)',[ones(1,size(yy,2))',aggyy(1,:)',yy(ind,:)']);
                        rrstat(ind,vss)=rtemp(2);
                        betadcdyalt(ind,vss)=betadcdytemp(2);
                        RR2_dcdy(ind,vss)=statstemp(1);
                        clear rtemp betadcdytemp statstemp
			[betadcdytemp,temp1,rtemp,temp2,statstemp]=regress(cc_cond(ind,:)',[ones(1,size(yy,2))',yy_cond(ind,:)']);
                        rrstat(ind,vss)=rtemp(2);
                        betadcdyalt1(ind,vss)=betadcdytemp(2);
                        RR2_dcdy_cond(ind,vss)=statstemp(1);
                        clear rtemp betadcdytemp statstemp

                        % for those with income gains versus
                        %those with income losses
                        YY=cc(ind,:)';
                        yy_hc(ind,:)=zeros(1,size(cc,2));
                        yy_hc(ind,rrhc)=yy(ind,rrhc);
                        oneone=ones(1,size(cc,2));
                        ones_hc(ind,:)=zeros(1,size(cc,2));
                        ones_hc(ind,rrhc)=oneone(rrhc);
                        XX=[ones(1,size(cc,2))',ones_hc(ind,:)',yy(ind,:)',yy_hc(ind,:)'];
                        beta_yhigh(:,ind,vss)=regress(YY,XX);
                        betadcdy_low(ind,vss)=beta_yhigh(3,ind,vss);
                        betadcdy_high(ind,vss)=beta_yhigh(4,ind,vss)+beta_yhigh(3,ind,vss);
                        
                        %2) conditional on village income
                        covtemp=[];
                        X=[ones(size(yy_cond(ind,:),2),1),yy_cond(ind,:)'];%
                        YY=cc_cond(ind,:)';
                        beta_cond(:,ind,vss)=inv(X'*X)*(X'*YY);
                        X=[ones(size(aggyy(1,:),2),1),aggyy'];%
                        YY=cc(ind,:)';
                        beta_agg(:,ind,vss)=inv(X'*X)*(X'*YY);
                        %for those with income gains versus
                        %those with income losses
                        YY_cond=cc_cond(ind,:)';
                        yy_hc_cond(ind,:)=zeros(1,size(cc_cond,2));
                        yy_hc_cond(ind,rrhc_cond)=yy_cond(ind,rrhc_cond);
                        oneone=ones(1,size(cc_cond,2));
                        ones_hc_cond(ind,:)=zeros(1,size(cc_cond,2));
                        ones_hc_cond(ind,rrhc_cond)=oneone(rrhc_cond);
                        XX_cond=[ones(1,size(cc_cond,2))',ones_hc_cond(ind,:)',yy_cond(ind,:)',yy_hc_cond(ind,:)'];
                        beta_yhigh_cond(:,ind,vss)=regress(YY_cond,XX_cond);

                        %3) conditional on group income
                        covtemp=[];
                        X=[ones(size(yy_groupcond(ind,:),2),1),yy_groupcond(ind,:)'];%
                        YY=cc_groupcond(ind,:)';
                        beta_groupcond(:,ind,vss)=inv(X'*X)*(X'*YY);
                        X=[ones(size(aggvillyy(ind,:),2),1),aggvillyy(ind,:)'];%
                        YY=cc(ind,:)';
                        beta_groupagg(:,ind,vss)=inv(X'*X)*(X'*YY);
                        %risk-sharing for those with income gains versus
                        %those with income losses
                        YY_groupcond=cc_groupcond(ind,:)';
                        yy_hc_groupcond(ind,:)=zeros(1,size(cc_groupcond,2));
                        yy_hc_groupcond(ind,rrhc_groupcond)=yy_groupcond(ind,rrhc_groupcond);
                        oneone=ones(1,size(cc_groupcond,2));
                        ones_hc_groupcond(ind,:)=zeros(1,size(cc_groupcond,2));
                        ones_hc_groupcond(ind,rrhc_groupcond)=oneone(rrhc_groupcond);
                        XX_groupcond=[ones(1,size(cc_groupcond,2))',ones_hc_groupcond(ind,:)',yy_groupcond(ind,:)',yy_hc_groupcond(ind,:)'];
                        beta_yhigh_groupcond(:,ind,vss)=regress(YY_groupcond,XX_groupcond);
                    end
                    var_dyagglog(1,vss)=var(aggyy);
                    R2_dcdy=mean(RR2_dcdy_cond(1:size(csim,1),vss));
                    
                    % =====================
                    % Matrix of moments
                    % =====================                    
                    % col 1-5 are parameters
                    Momentsest1(1:5,ic)=[vss,beta,rho,rho1vec(incproc),varnu(incproc)]';
                    % 6-11 relative variance dc/dy, dc(dy>0)/dc(dy<0);  unconditional (6:7), conditional on village income
                    % (8:9), conditonal on group income (10:11)
                    Momentsest1(6:11,ic)=[mean(vardclog(1:size(csim,1),vss))/mean(vardylog(1:size(csim,1),vss)), mean(reldcvar(1:size(csim,1),vss))...
                    mean(vardclog_cond(1:size(csim,1),vss))/mean(vardylog_cond(1:size(csim,1),vss)), mean(reldcvar_cond(1:size(csim,1),vss))...
                    mean(vardclog_groupcond(1:size(csim,1),vss))/mean(vardylog_groupcond(1:size(csim,1),vss)), mean(reldcvar_groupcond(1:size(csim,1),vss))]';
                    % 12-14: beta, beta high, beta low unconditional
                    Momentsest1(12:14,ic)=[mean(betadcdy(1:size(csim,1),vss)),mean(betadcdy_high(1:size(csim,1),vss)),mean(betadcdy_low(1:size(csim,1),vss))]';
                    % 13_16: beta, beta high, beta low conditional on
                    % village income
                    Momentsest1(15:17,ic)=[mean(beta_cond(2,1:size(csim,1),vss),2),mean(beta_yhigh_cond(3,1:size(csim,1),vss),2)+mean(beta_yhigh_cond(4,1:size(csim,1),vss),2),mean(beta_yhigh_cond(3,1:size(csim,1),vss),2)]';
                    % 18:20: beta, beta high, beta low conditional on
                    % group income
                    Momentsest1(18:20,ic)=[mean(beta_groupcond(2,1:size(csim,1),vss),2),mean(beta_yhigh_groupcond(3,1:size(csim,1),vss),2)+mean(beta_yhigh_groupcond(4,1:size(csim,1),vss),2),mean(beta_yhigh_groupcond(3,1:size(csim,1),vss),2)]';
                    %21:22 are variances of dclog, dylog
                    Momentsest1(21:22,ic)=[mean(vardclog(1:size(csim,1),vss)),mean(vardylog(1:size(csim,1),vss))];
                    %23:28, relative variance dc/dy (dy>0) and var dc/dy(dy<0);
                    %unconditinoal 23:24, conditional on village income 25:26,
                    %conditional on group income 27:28
                    Momentsest1(23:24,ic)=[mean(vardc_dypos(1:size(csim,1),vss));mean( vardc_dyneg(1:size(csim,1),vss))];
                    Momentsest1(25:26,ic)=[mean(vardc_dypos_cond(1:size(csim,1),vss));mean(  vardc_dyneg_cond(1:size(csim,1),vss))];
                    Momentsest1(27:28,ic)=[mean(vardc_dypos_groupcond(1:size(csim,1),vss));mean( vardc_dyneg_groupcond(1:size(csim,1),vss))];
                    Momentsest1(30,ic)=[mean( vardylog_cond(1:size(csim,1),vss))];
                    Momentsest1(31,ic)=mean(RR2_dcdy(1:size(csim,1),vss));
                    Momentsest1(32,ic)=mean(vardclog(1:size(csim,1),vss));
                    Momentsest1(33,ic)=mean(vardclog_cond(1:size(csim,1),vss));
 		   Momentsest1(34,ic)=var_dc_t_fe(vss);
 		    Momentsest1(35,ic)=mean(RR2_dcdy_cond(1:size(csim,1),vss));
                    %%% and some other model-dependent measures
                    if coaldev<2
                        Momentsest1(29,ic)=CRITGAP(vss);
                    elseif coaldev==2
                        Momentsest1(29,ic)=phi;
                    elseif coaldev==3
                        Momentsest1(29,ic)=nan;
                    end
                end
                Momentsest(:,:,Qbeta,Qrho,vss,equal+1)=[Momentsest1];
                toc
            end

                cd(outputpath)
                eval(['save ', filetosaveest, int2str(equal) '.mat'])
                cd(programpath)
        end
    end
    end
 %%%%%%%%%%%%%%%%%%%%%%%%%%%do your own Momentsexplore down here  
 analytic=1
 outlier=0 
 relvarincscale=1;
 incproc=1
village=1;                         % village indicator: 1 - Aurepalle; 2-Kanzara; 3-Shirapur
coaldev=1;                         % set to 1 for coalitional deviations, to 2 for self-insurance (SI), to 0 for individual deviations (ID), to 3 for full-insurance with coalitional deviations
Tessa1Tobi2=1;                     % sets paths depending on user below
unconditional=0;                   % set this to 1 for unconditional income data, to 0 for income data conditional on time fixed effets / demeaned data period-by-period
fixedeffect=1;                     % set this to 1 for an income process with fixed effects
one_per_aut=0;                     % set this to 1 if you want deviations to imply one period of autarky
ytaxrate=0;                        % linear tax on income
scatter = 0;		           % set to 1 to just draw a scatter plot
Rov_sb = 0; 
 if Tessa1Tobi2==1
    outputpath='C:\Users\tbold\Dropbox\Tessa and Tobi\Documents\Latex';
        datapath='C:\Users\tbold\Dropbox\Tessa and Tobi\JEEA Replication\Data/Output';
        
else
    outputpath='C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\JEEA Replication\Documents\Latex';
    datapath='C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\JEEA Replication\Data\Output';
end   

        datafile=['icrisatmoments' int2str(outlier) ];
   
    
        cd(datapath)
  
    eval(['load ' datafile '.out'])
    eval(['datamoments=' datafile ';']);
    eval(['clear ' datafile])
    
    
    Nom=4;
    
    bootstrapvar=datamoments(1:Nom,4+(village-1)*Nom+1:4+(village)*Nom);
    analyticvar=datamoments(1:Nom,16+(village-1)*Nom+1:16+(village)*Nom);
    analyticmoments=datamoments(1:Nom,1:3);
    
    analyticrawmoments=datamoments(1:Nom,28+1:28+3);
    bootstraprawvar=datamoments(1:Nom,28+4+(village-1)*Nom+1:28+4+(village)*Nom);
    analyticrawvar=datamoments(1:Nom,28+16+(village-1)*Nom+1:28+16+(village)*Nom);
    %%%here you replace the exact standard errors (the variances are
    %%%bootstrapped in STATA, the coefficients are delta-method
    
    VCV2data=bootstrapvar;
    VCV2rawdata=bootstraprawvar;
    
    if analytic==1
        VCV2data(2,2)=analyticvar(2,2);
        VCV2data(4,4)=analyticvar(4,4);
        VCV2rawdata(2,2)=analyticrawvar(2,2);
        VCV2rawdata(4,4)=analyticrawvar(4,4);
    end
        
   
Momvec2data(1,1:4)=[0.30 0.21 -0.08 -0.41];
Vardylog_cond=[0.336;0.162;0.268];
        Momentssum=[];
   
          
    
        for equal=0:1
           deltaset=[find(isfinite(squeeze(Momentsest(1,1,:,1,1,equal+1))) &squeeze(Momentsest(1,1,:,1,1,equal+1))>0 )]';
         
            for inddelta=deltaset
                for indrho=1
                    for j=1
                        vss_set=[find((squeeze(Momentsest(1,j,inddelta,indrho,:,equal+1)))>0 & isfinite(squeeze(Momentsest(1,j,inddelta,indrho,:,equal+1))))]';
                        vss_max=max(vss_set)
                                if coaldev==0
                                vss_set=vss_max;
                                end
                        for vss=vss_max
                            Momentstable(1:3,inddelta,equal+1)=Momentsest(1:3,j,inddelta,indrho,vss,equal+1);
                            if relvarincscale==1
                            Momentstable(4:7,inddelta,equal+1)=[Momentsest([8 15],j,inddelta,indrho,vss,equal+1) ; Momentsest(9,j,inddelta,indrho,vss,equal+1)/Vardylog_cond(village) ;(Momentsest([16],j,inddelta,indrho,vss,equal+1)-Momentsest([17],j,inddelta,indrho,vss,equal+1))];
                            end
                            Momentstable(8,inddelta,equal+1)=(Momentsest(21,j,inddelta,indrho,vss,equal+1)-Momentsest(21,1,inddelta,indrho,vss,equal+1))./Momentsest(21,j,inddelta,indrho,vss,equal+1);%2*maxcerr*(j-1)/(size(moments,2)-1)/moments(28,j);
                            Momentstable(9,inddelta,equal+1)=2*Momentsest(5,j,inddelta,indrho,vss,equal+1);%Momentsest(22,j,inddelta,indrho,vss);
                            Momentstable(1,inddelta,equal+1)=Momentstable(1,inddelta,equal+1)+1;
                            
                             
                                
                            momentsdev=Momentstable(4:7,inddelta,equal+1)- Momvec2data(1,1:4)';
                            %%%2,4, d:diagonal, 
                            
                            
                            
                            W4(inddelta,equal+1)=momentsdev'*inv(VCV2data(1:4,1:4))*momentsdev;
                            W4d(inddelta,equal+1)=momentsdev'*inv(diag(diag(VCV2data(1:4,1:4))))*momentsdev;
                            
                               
                            
                        end
                        
                           
                    end
                   
                        
               end
            end
        end
        
%%Estimation (for unequal sharing)
  WW=W4d(:,1);
  [CC,II]=min(WW);
  solution_4(:,1)=[Momentstable(1:9,II,1); CC];
     
 solution_4(:,1)
%%Estimation (for equal sharing)         
 WW=W4d(:,2);
  [CC,II]=min(WW);
  solution_4(:,2)=[Momentstable(1:9,II,2); CC];
  
 solution_4(:,2)